rm(list=ls(all=TRUE))

### data file created by running script_01_relicate_original_results.R
load("Data/Temp Files/BIW_replicated_datasets.RData")

# ### add fraction of copartisans deviating
# da <- da %>% 
#   mutate(. , roll_call_party = paste0(roll_call , r_party)) %>% 
#   .[order(.$roll_call_party,  .$bioname),]
# 
# da <- da[, `:=` (oparty_fracdeviate_loo = (sum(na.omit(deviate)) - deviate)/(length(na.omit(deviate))-1)) , 
#          by = roll_call_party ]
# da$oparty_fracdeviate_loo_M <- as.numeric(is.na(da$oparty_fracdeviate_loo))
# da$oparty_fracdeviate_loo[is.na(da$oparty_fracdeviate_loo)] <- 1


### roll call fixed effects
col1 <- feols(deviate ~ rank_p  | icpsr + roll_call, data = subset(da2 , congress < 113 &!cast_code%in%c(2,3,4,5,9) & r_party!="other"))
col2 <- feols(deviate ~ rank_p  | interaction(icpsr,congress) + roll_call, data = subset(da2 , congress < 113 &!cast_code%in%c(2,3,4,5,9) & r_party!="other")) 

col3 <- feols(deviate ~ rank_p_overall  | icpsr + roll_call, data = subset(da2 , congress < 113 &!cast_code%in%c(2,3,4,5,9) & r_party!="other")) 
col4 <- feols(deviate ~ rank_p_overall  | interaction(icpsr,congress) + roll_call, data = subset(da2 , congress < 113 &!cast_code%in%c(2,3,4,5,9) & r_party!="other")) 

tableb1a <- list(col1 , col2 , col3 , col4 )

### create txt file with results
sink("Generated Files/Replicated Table B1a.txt")
column <- 1
for(i in tableb1a){
  cat("Column ",column,"\n")
  summary(i , cluster = "congress") %>% 
    print
  cat("R2 = " , round(r2(i)[2],3),"\n\n")
  column <- column + 1
}
sink(file = NULL)


### roll call by party fixed effects
col1 <- feols(deviate ~ rank_p  | icpsr + roll_call_party, data = subset(da2 , congress < 113 &!cast_code%in%c(2,3,4,5,9) & r_party!="other"))
col2 <- feols(deviate ~ rank_p  | interaction(icpsr,congress) + roll_call_party, data = subset(da2 , congress < 113 &!cast_code%in%c(2,3,4,5,9) & r_party!="other")) 

col3 <- feols(deviate ~ rank_p_overall  | icpsr + roll_call_party, data = subset(da2 , congress < 113 &!cast_code%in%c(2,3,4,5,9) & r_party!="other")) 
col4 <- feols(deviate ~ rank_p_overall  | interaction(icpsr,congress) + roll_call_party, data = subset(da2 , congress < 113 &!cast_code%in%c(2,3,4,5,9) & r_party!="other")) 

tableb1b <- list(col1 , col2 , col3 , col4 )

### create txt file with results
sink("Generated Files/Replicated Table B1b.txt")
column <- 1
for(i in tableb1b){
  cat("Column ",column,"\n")
  summary(i , cluster = "congress") %>% 
    print
  cat("R2 = " , round(r2(i)[2],3),"\n\n")
  column <- column + 1
}
sink(file = NULL)


#####################################################################
#############    Randomization inference - table 1   ##############
#####################################################################
### results from these simulations also used in Table 5 and associated figures
nsims <- 10000
output3 <- output1 <- output2 <- data.frame(matrix(NA , nrow = nsims , ncol = 8))
output3rc <- output1rc <- output2rc <- data.frame(matrix(NA , nrow = nsims , ncol = 8))
output1_wrongri <- data.frame(matrix(NA , nrow = nsims , ncol = 4))
### j: 1 is for original timeframe; 2 is for new data; 3 is for all data
### k: 1/2 is for for roll-call specific vote order or overall vote order; 
### k: 3 for wrong ri - repeating error in original BIW code of dropping non-Ds and Rs
  ### before permuting names

### have to use full data to permute names, to directly replicate original analysis
names.unique <- unique(da$bioname)

set.seed(12345)
for(i in 1:nsims){
  temp <- data.frame(bioname = names.unique,
                     name2 = sample(names.unique , length(names.unique) , replace = FALSE)) %>% as.data.table()
  da_sim <- join(da , temp , by = "bioname")

  for(k in 1:3){
    if(k==1){
      da_sim <- da_sim[order(da_sim$roll_call , da_sim$name2),]
      da_sim2 <- da_sim[, `:=` (rank_p_sim = (seq_along(name2)-1)/(length(name2)-1)) , by = roll_call ]
    }
    if(k==2){
      da_sim2 <- subset(da_sim , !cast_code%in%c(9))
      da_sim2 <- da_sim2[order(da_sim2$roll_call , da_sim2$name2),]
      da_sim2 <- da_sim2[, `:=` (rank_p_sim = (seq_along(name2)-1)/(length(name2)-1)) , by = roll_call ]
    }
    
    if(k==3){
      da_sim2 <- subset(da_sim , !cast_code%in%c(2,3,4,5,9) & r_party!="other" & !is.na(deviate))
      da_sim2 <- da_sim2[order(da_sim2$roll_call , da_sim2$name2),]
      da_sim2 <- da_sim2[, `:=` (rank_p_sim = (seq_along(name2)-1)/(length(name2)-1)) , by = roll_call ]
    }
    
    #much smaller standard errors
    for(j in 1:3){
      if(j==1) da_sim3 <- subset(da_sim2 , congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other" & !is.na(deviate))
      if(j==2 & k!=3) da_sim3 <- subset(da_sim2 , congress >= 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other" & !is.na(deviate))
      if(j==3 & k!=3) da_sim3 <- subset(da_sim2 , !cast_code%in%c(2,3,4,5,9) & r_party!="other" & !is.na(deviate))
      
      tttemp <- matrix(NA , nrow = 1 , ncol = 8)
      ttemp <- feols(deviate ~ rank_p_sim  | name2, data = da_sim3)
      ttemp1 <- summary(ttemp , cluster = "congress") %>% coeftable(.)
      tttemp[1,1:4] <- unlist(ttemp1)
      if(k!=3){
        ttemp <- feols(deviate ~ rank_p_sim  | interaction(name2 , congress), data = da_sim3)
        ttemp1 <- summary(ttemp , cluster = "congress") %>% coeftable(.)
        tttemp[1,5:8] <- unlist(ttemp1)
      }  

      if(j==1 & k==1) output1[i,] <- tttemp
      if(j==2 & k==1) output2[i,] <- tttemp
      if(j==3 & k==1) output3[i,] <- tttemp
      if(j==1 & k==2) output1rc[i,] <- tttemp
      if(j==2 & k==2) output2rc[i,] <- tttemp
      if(j==3 & k==2) output3rc[i,] <- tttemp
      if(j==1 & k==3) output1_wrongri[i,] <- tttemp[,1:4]
    }
  }
  cat("\r", i, "of", nsims) ### How far loop has progressed
}


save(output1,output2,output3,output1rc,output2rc,output3rc,output1_wrongri,
     file = "Data/Temp Files/BIW_replicated_RI_sims.RData")

################################################################################
####################### LOAD DATA AND START FROM HERE ##########################
################################################################################
load(file = "Data/Temp Files/BIW_replicated_RI_sims.RData")


################################################################################
########################### For Plots - skip if you want #######################
################################################################################

if(FALSE){
library(myriad)
library(extrafont)
require(ggokabeito)
loadfonts(device = "win")
theme_figs <- function(){
  theme_myriad_semi(base_family = fonts()[1],
                    subtitle_family = fonts()[1],
                    caption_family = fonts()[1]) +
    theme(
      plot.background = element_rect(color = "white"),
      plot.title = element_text(size = rel(1.5)),
      plot.subtitle = element_text(size = rel(1.2)),
      axis.title.x = element_text(size = rel(1.25)),
      axis.title.y = element_text(size = rel(1.25)),
      axis.text.y.left = element_text(size = rel(1.2)),
      axis.text.y.right = element_text(size = rel(0.9)))
}
theme_set(theme_figs())
## Colors for manual use as needed
my_oka <- palette_okabe_ito(order = c(1, 2, 3, 5, 7),
                            alpha = NULL)
my_oka2 <- palette_okabe_ito(order = c(1, 2, 3, 5, 7),
                             alpha = rep(0.25))
}
################################################################################
################################################################################
pdf(
  "Generated Files/Figure_5.pdf",
  height = 5,
  width = 8,
  onefile = FALSE
)
titles <- c("Leg FE","Leg x Session FE")
par(mfcol = c(1,1))
ylim2 <- 4000
hist(output1rc[,4] , breaks = 20 ,
     ylim = c(0,ylim2),
     xlab = "Est. p-value", ylab = "Frequency",
     main = "")
dev.off()
mean(output1rc[,4] <= 0.05)
mean(output1rc[,4] <= 0.01)
mean(output1rc[,4] <= 0.001)

### Table 1 is created via these simulations
table1 <- matrix(NA , nrow = 3 , ncol = 4)


#############################################################
### reproduce CDFs from simulations. Ranking within each roll call, like BIW
### also Table 1
### original data
estimate <- feols(deviate ~ rank_p  | icpsr, data = subset(da2 , congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
sd(output1rc[,1])
mean(output1rc[,1] <= estimate[1,1] )
mean(abs(output1rc[,1]) >= abs(estimate[1,1]) )
mean(output1_wrongri[,1] <= estimate[1,1] )

### can use deltax variable to interpret estimated effects in terms of # of spots in vote order
#deltax <- 100*1/99
deltax <- 1

g1 <- ggplot(output1rc, aes(X1*deltax)) +
  stat_ecdf(geom = "step")+
  labs(#title="Leg FE",
    subtitle = "Original Data",
    y = "Cumulative Probability", x="Alphabetical Rank Estimate")+
  #3theme_classic()+
  geom_hline(yintercept = mean(output1rc$X1 <=  -0.17023),
             col = "black" , linetype="longdash") +
  geom_text(data = data.frame(x1 =  -0.48,
                              y1 = mean(output1rc$X1 <=  -0.17023 + 0.03),
                              label1 = round(mean(output1rc$X1 <=  -0.17023),3)),
            aes(x = x1, y = y1, label = label1))+
  scale_y_continuous(limits = c(0,1),breaks = seq(0,1,0.1))+
  scale_x_continuous(limits = c(-0.5,0.5),breaks = seq(-0.5,0.5,0.25))+
  geom_vline(xintercept = -0.17023*deltax,
             col = "black" , linetype="longdash") #+

pdf(
  "Generated Files/Figure_4.pdf",
  height = 4,
  width = 5,
  onefile = FALSE
)
g1 + labs(subtitle = NULL)
dev.off()





estimate <-feols(deviate ~ rank_p  | interaction(icpsr,congress) , data = subset(da2 , congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
sd(output1rc[,5])
mean(output1rc[,5] <= estimate[1,1] )
mean(abs(output1rc[,5]) >= abs(estimate[1,1]))

g2 <- ggplot(output1rc, aes(X5*deltax)) +
  stat_ecdf(geom = "step")+
  labs(#title="Leg by Congress FE",
    subtitle = "Original Data",
    y = "Cumulative Probability", x="Alphabetical Rank Estimate")+
  # theme_classic()+
  geom_hline(yintercept = mean(output1rc$X5 <=  estimate[1,1]),
             col = "black" , linetype="longdash") +
  geom_text(data = data.frame(x1 =  -.48,
                              y1 = mean(output1rc$X5 <=  estimate[1,1] + 0.02),
                              label1 = round(mean(output1rc$X5 <=  estimate[1,1]),3)),
            aes(x = x1, y = y1, label = label1))+
  scale_y_continuous(limits = c(0,1),breaks = seq(0,1,0.1))+
  scale_x_continuous(limits = c(-0.5,0.5),breaks = seq(-0.5,0.5,0.25))+
  geom_vline(xintercept = estimate[1,1]*deltax,
             col = "black" , linetype="longdash")


### new data
estimate <-feols(deviate ~ rank_p  | icpsr, data = subset(da2 , congress >= 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
sd(output2rc[,1])
mean(output2rc[,1] <= estimate[1,1] )
mean(abs(output2rc[,1]) >= abs(estimate[1,1]))


g1new <- ggplot(output2rc, aes(X1*deltax)) +
  stat_ecdf(geom = "step")+
  labs(#title="",
    subtitle = "New Data",
    y = "", x="Alphabetical Rank Estimate")+
  # theme_classic()+
  geom_hline(yintercept = mean(output2rc$X1 <=  estimate[1,1]),
             col = "black" , linetype="longdash") +
  geom_text(data = data.frame(x1 =  -0.48,
                              y1 = mean(output2rc$X1 <=  estimate[1,1] + 0.015),
                              label1 = round(mean(output2rc$X1 <=  estimate[1,1]),3)),
            aes(x = x1, y = y1, label = label1))+
  scale_y_continuous(limits = c(0,1),breaks = seq(0,1,0.1))+
  scale_x_continuous(limits = c(-0.5,0.5),breaks = seq(-0.5,0.5,0.25))+
  geom_vline(xintercept = estimate[1,1]*deltax,
             col = "black" , linetype="longdash")



estimate <-feols(deviate ~ rank_p  | interaction(icpsr,congress) , data = subset(da2 , congress >= 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
sd(output2rc[,5])
mean(output2rc[,5] <= estimate[1,1] )
mean(abs(output2rc[,5]) >= abs(estimate[1,1]))


### all data
estimate <-feols(deviate ~ rank_p  | icpsr, data = subset(da2 , !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
sd(output3rc[,1])
mean(output3rc[,1] <= estimate[1,1] )
mean(abs(output3rc[,1]) >= abs(estimate[1,1]) )

g1all <- ggplot(output3rc, aes(X1*deltax)) +
  stat_ecdf(geom = "step")+
  labs(#title="",
    subtitle = "All Data",
    y = "", x="Alphabetical Rank Estimate")+
  #  theme_classic()+
  geom_hline(yintercept = mean(output3rc$X1 <=  estimate[1,1]),
             col = "black" , linetype="longdash") +
  geom_text(data = data.frame(x1 =  -0.48,
                              y1 = mean(output3rc$X1 <=  estimate[1,1] + 0.025),
                              label1 = format(round(mean(output3rc$X1 <=  estimate[1,1]),3),nsmall=3)),
            aes(x = x1, y = y1, label = label1))+
  scale_y_continuous(limits = c(0,1),breaks = seq(0,1,0.1))+
  scale_x_continuous(limits = c(-0.5,0.5),breaks = seq(-0.5,0.5,0.25))+
  geom_vline(xintercept = estimate[1,1]*deltax,
             col = "black" , linetype="longdash")

estimate <-feols(deviate ~ rank_p  | interaction(icpsr,congress) , data = subset(da2 ,  !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
sd(output3rc[,5])
mean(output3rc[,5] <= estimate[1,1] )
mean(abs(output3rc[,5]) >= abs(estimate[1,1]) )

### now create columns 3 and 4 of Table 1
estimate <- feols(deviate ~ rank_p_overall  | icpsr, data = subset(da2 , congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(abs(output1[,1]) >= abs(estimate[1,1]) )

estimate <- feols(deviate ~ rank_p_overall  | icpsr, data = subset(da2 , congress >= 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(abs(output2[,1]) >= abs(estimate[1,1]) )

estimate <- feols(deviate ~ rank_p_overall  | icpsr, data = subset(da2 ,  !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(abs(output3[,1]) >= abs(estimate[1,1]) )

# column 4
estimate <- feols(deviate ~ rank_p_overall  | interaction(icpsr,congress), data = subset(da2 , congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(abs(output1[,5]) >= abs(estimate[1,1]) )

estimate <- feols(deviate ~ rank_p_overall  | interaction(icpsr,congress), data = subset(da2 , congress >= 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(abs(output2[,5]) >= abs(estimate[1,1]) )

estimate <- feols(deviate ~ rank_p_overall  | interaction(icpsr,congress), data = subset(da2 ,  !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(abs(output3[,5]) >= abs(estimate[1,1]) )


require(ggpubr)
pdf(
  "Generated Files/Figure_2.pdf",
  height = 4,
  width = 12,
  onefile = FALSE
)
ggarrange(g1, g1new , g1all,
          nrow = 1, ncol = 3)
dev.off()




################################################################################
#############################    Create Figure 1   #############################
################################################################################

nsims <- 1000
ri.windows <- matrix(NA , nrow = nsims , ncol = 63)

names.unique <- unique(da$bioname)

set.seed(12345)
for(i in 1:nsims){
  temp <- data.frame(bioname = names.unique,
                     name2 = sample(names.unique , length(names.unique) , replace = FALSE)) %>% as.data.table()
  datemp <- join(da , temp , by = "bioname")
  
  datempa <- subset(datemp , !cast_code%in%c(9))
  datempa <- datempa[order(datempa$roll_call , datempa$name2),]
  datempa <- datempa[, `:=` (rank_p_rc = (seq_along(name2)-1)/(length(name2)-1)) , by = roll_call ]
  
  datempb <- subset(datempa , !cast_code%in%c(2,3,4,5,9) & r_party!="other" & !is.na(deviate))
  for(k in seq(35,97 , 1)){
    tempdata <- subset(datempb , congress >= k & congress < k + 20)
    #  regfe2 <- feols(deviate ~ rank_p  | interaction(icpsr,congress) ,  data = tempdata)
    ri.windows[i,match(k,seq(35,112 , 1))] <- feols(deviate ~ rank_p_rc  | name2 ,  data = tempdata) %>% coeftable %>% .[1,1]
  }
  cat("\r", i, "of", nsims) ### How far loop has progressed
}


save(ri.windows ,
     file = "Data/Temp Files/us_sims_windows.RData")

load("Data/Temp Files/us_sims_windows.RData")
dim(ri.windows)

temp <- data.frame(matrix(NA , nrow = 63 , ncol = 2))
tempdata2 <- da2 %>%
  subset(. , !cast_code%in%c(2,3,4,5,9) & !is.na(deviate) & r_party!="other")

for(i in seq(35,97 , 1)){
  tempdata <- subset(tempdata2 , congress >= i & congress < i + 20)
  regfe2 <- feols(deviate ~ rank_p  | icpsr, data =tempdata)
  temp[match(i,seq(35,112 , 1)),1:2] <- summary(regfe2, cluster = "congress")[c('coefficients','se')] %>% unlist()
  
}

temp$se_ri <- apply(ri.windows , 2 , sd)

temp$lower <- temp$X1 - 1.65*temp$se_ri
temp$upper <- temp$X1 + 1.65*temp$se_ri
temp$cong <- c(55:117)

pdf(
  "Generated Files/Figure_1.pdf",
  height = 4,
  width = 10,
  onefile = FALSE
)
ggplot(data = subset(temp , cong <= 112)) +
  geom_pointrange(mapping = aes(x = cong , y = X1 , ymin = lower, ymax = upper))+
  geom_pointrange(data = subset(temp , cong > 112) , mapping = aes(x = cong , y = X1 , ymin = lower, ymax = upper),
                  col = "red",shape = 15)+
  xlab("Congress (end of window)")+ ylab("Alphabetical Rank Estimate")
dev.off()



#####################################################################
#############    Randomization inference - close v lopsided   ##############
#####################################################################
nsims <- 10000
output_all_rc <- output_close_rc <- output_lopsided_rc <- output_all <- output_close <- 
  output_lopsided <- data.frame(matrix(NA , nrow = nsims , ncol = 12))

### j: 1 is for original timeframe; 2 is for new data; 3 is for all data
### k: 1/2 is for for roll-call specific vote order or overall vote order; 
### k: 3 for wrongri repeating error in original BIW code of dropping non-Ds and Rs
### before permuting names

### have to use full data to permute names, to directly replicate original analysis
da_temp1 <- subset(da , congress < 113) 
# da_temp1 <- da_temp1[, `:=` (rank_p_sim = (seq_along(name2)-1)/(length(name2)-1),
#   ) , by = roll_call ]
names.unique <- unique(da_temp1$bioname)

set.seed(12345)
for(i in 1:nsims){
  temp <- data.frame(bioname = names.unique,
                     name2 = sample(names.unique , length(names.unique) , replace = FALSE)) %>% as.data.table()
  da_sim <- join(da_temp1 , temp , by = "bioname")
  
  for(k in 1:2){
    if(k==1){
      da_sim <- da_sim[order(da_sim$roll_call , da_sim$name2),]
      da_sim2 <- da_sim[, `:=` (rank_p_sim = (seq_along(name2)-1)/(length(name2)-1)) , by = roll_call ]
    }
    if(k==2){
      da_sim2 <- subset(da_sim , !cast_code%in%c(9))
      da_sim2 <- da_sim2[order(da_sim2$roll_call , da_sim2$name2),]
      da_sim2 <- da_sim2[, `:=` (rank_p_sim = (seq_along(name2)-1)/(length(name2)-1)) , by = roll_call ]
    }
    
    # all, close, lopsided
    for(j in 1:3){
      if(j==1) da_sim3 <- subset(da_sim2 , !cast_code%in%c(2,3,4,5,9) & r_party!="other" & !is.na(deviate))
      if(j==2) da_sim3 <- subset(da_sim2 , r_close==0 & !cast_code%in%c(2,3,4,5,9) & r_party!="other" & !is.na(deviate))
      if(j==3) da_sim3 <- subset(da_sim2 , r_close==1 & !cast_code%in%c(2,3,4,5,9) & r_party!="other" & !is.na(deviate))
      
      tttemp <- matrix(NA , nrow = 1 , ncol = 12)
      ttemp <- feols(deviate ~ rank_p_sim  | name2, data = da_sim3)
      ttemp1 <- summary(ttemp , cluster = "congress") %>% coeftable(.)
      tttemp[1,1:4] <- unlist(ttemp1)

      ttemp <- feols(deviate ~ rank_p_sim  | interaction(name2,congress), data = da_sim3)
      ttemp1 <- summary(ttemp , cluster = "congress") %>% coeftable(.)
      tttemp[1,5:8] <- unlist(ttemp1)
      
      ttemp <- feols(deviate ~ rank_p_sim + oparty_fracdeviate_loo + oparty_fracdeviate_loo_M | interaction(name2,congress), data = da_sim3)
      ttemp1 <- summary(ttemp , cluster = "congress") %>% coeftable(.)
      tttemp[1,9:12] <- unlist(ttemp1[1,])
      
      if(j==1 & k==1) output_all[i,] <- tttemp
      if(j==2 & k==1) output_lopsided[i,] <- tttemp
      if(j==3 & k==1) output_close[i,] <- tttemp
      if(j==1 & k==2) output_all_rc[i,] <- tttemp
      if(j==2 & k==2) output_lopsided_rc[i,] <- tttemp
      if(j==3 & k==2) output_close_rc[i,] <- tttemp
    }
  }
    cat("\r", i, "of", nsims) ### How far loop has progressed
}


save(output_all , output_close , output_lopsided,output_all_rc , output_close_rc , output_lopsided_rc,
     file = "Data/Temp Files/BIW_replicated_RI_sims_close_lopsided.RData")


### Create Table 7 - note slight discrepancies with Table 1, as RI relies on distinct
### simulations
load("Data/Temp Files/BIW_replicated_RI_sims_close_lopsided.RData")

da2 <- subset(da , !cast_code%in%c(9))
da2 <- da2[order(da2$roll_call , da2$bioname),]
da2 <- da2[, `:=` (rank_p = (seq_along(bioname)-1)/(length(bioname)-1)) , by = roll_call ]


### all votes - 
estimate <- feols(deviate ~ rank_p  | icpsr, data = subset(da2 , congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(output_all_rc[,1] <= estimate[1,1] )
mean(abs(output_all_rc[,1]) >= abs(estimate[1,1]) )
### close - off due to sampling variability

estimate <- feols(deviate ~ rank_p  | interaction(icpsr,congress) ,  data = subset(da2 , congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(output_all_rc[,5] <= estimate[1,1] )
mean(abs(output_all_rc[,5]) >= abs(estimate[1,1]) )

estimate <- feols(deviate ~ rank_p + oparty_fracdeviate_loo + oparty_fracdeviate_loo_M | interaction(icpsr,congress) ,  data = subset(da2 , congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(output_all_rc[,9] <= estimate[1,1] )
mean(abs(output_all_rc[,9]) >= abs(estimate[1,1]) )


estimate <- feols(deviate ~ rank_p_overall  | icpsr, data = subset(da2 , congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(output_all[,1] <= estimate[1,1] )
mean(abs(output_all[,1]) >= abs(estimate[1,1]) )


estimate <- feols(deviate ~ rank_p_overall  | interaction(icpsr,congress), data = subset(da2 , congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(output_all[,5] <= estimate[1,1] )
mean(abs(output_all[,5]) >= abs(estimate[1,1]) )


estimate <- feols(deviate ~ rank_p_overall + oparty_fracdeviate_loo + oparty_fracdeviate_loo_M | interaction(icpsr,congress) ,  data = subset(da2 , congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(output_all[,9] <= estimate[1,1] )
mean(abs(output_all[,9]) >= abs(estimate[1,1]) )


### close votes
estimate <- feols(deviate ~ rank_p  | icpsr, data = subset(da2 , r_close==1 & congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(output_close_rc[,1] <= estimate[1,1] )
mean(abs(output_close_rc[,1]) >= abs(estimate[1,1]) )

estimate <- feols(deviate ~ rank_p  | interaction(icpsr,congress) ,  data = subset(da2 ,r_close==1 &  congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(output_close_rc[,5] <= estimate[1,1] )
mean(abs(output_close_rc[,5]) >= abs(estimate[1,1]) )

estimate <- feols(deviate ~ rank_p + oparty_fracdeviate_loo + oparty_fracdeviate_loo_M | interaction(icpsr,congress) ,  data = subset(da2 ,r_close==1 &  congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(output_close_rc[,9] <= estimate[1,1] )
mean(abs(output_close_rc[,9]) >= abs(estimate[1,1]) )



estimate <- feols(deviate ~ rank_p_overall  | icpsr, data = subset(da2 , r_close==1 & congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(output_close[,1] <= estimate[1,1] )
mean(abs(output_close[,1]) >= abs(estimate[1,1]) )


estimate <- feols(deviate ~ rank_p_overall  | interaction(icpsr,congress), data = subset(da2 , r_close==1 & congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(output_close[,5] <= estimate[1,1] )
mean(abs(output_close[,5]) >= abs(estimate[1,1]) )


estimate <- feols(deviate ~ rank_p_overall + oparty_fracdeviate_loo + oparty_fracdeviate_loo_M | interaction(icpsr,congress) ,  data = subset(da2 ,r_close==1 &  congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(output_close[,9] <= estimate[1,1] )
mean(abs(output_close[,9]) >= abs(estimate[1,1]) )


### lopsided votes
estimate <- feols(deviate ~ rank_p  | icpsr, data = subset(da2 , r_close==0 & congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(output_lopsided_rc[,1] <= estimate[1,1] )
mean(abs(output_lopsided_rc[,1]) >= abs(estimate[1,1]) )

estimate <- feols(deviate ~ rank_p  | interaction(icpsr,congress) ,  data = subset(da2 ,r_close==0 &  congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(output_lopsided_rc[,5] <= estimate[1,1] )
mean(abs(output_lopsided_rc[,5]) >= abs(estimate[1,1]) )

estimate <- feols(deviate ~ rank_p + oparty_fracdeviate_loo + oparty_fracdeviate_loo_M | interaction(icpsr,congress) ,  data = subset(da2 ,r_close==0 &  congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(output_lopsided_rc[,9] <= estimate[1,1] )
mean(abs(output_lopsided_rc[,9]) >= abs(estimate[1,1]) )


estimate <- feols(deviate ~ rank_p_overall  | icpsr, data = subset(da2 , r_close==0 & congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(output_lopsided[,1] <= estimate[1,1] )
mean(abs(output_lopsided[,1]) >= abs(estimate[1,1]) )


estimate <- feols(deviate ~ rank_p_overall  | interaction(icpsr,congress), data = subset(da2 , r_close==0 & congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(output_lopsided[,5] <= estimate[1,1] )
mean(abs(output_lopsided[,5]) >= abs(estimate[1,1]) )

estimate <- feols(deviate ~ rank_p_overall + oparty_fracdeviate_loo + oparty_fracdeviate_loo_M | interaction(icpsr,congress) ,  data = subset(da2 ,r_close==0 &  congress < 113 & !cast_code%in%c(2,3,4,5,9) & r_party!="other")) %>%
  summary(. , cluster = "congress") %>% coeftable()
estimate
mean(output_lopsided[,9] <= estimate[1,1] )
mean(abs(output_lopsided[,9]) >= abs(estimate[1,1]) )




